#ifndef _H_EFLOAT
#define _H_EFLOAT
#include "heads.hpp"

namespace chen {

class EFloat;

inline bool isNan (EFloat a) ;

class EFloat {
public:
    Float value, l, h;
    inline EFloat() {
        value = l = h = 0;
    }
    inline explicit EFloat (const Float & v): value(v), l(v), h(v) {}
    inline EFloat (const Float & v, const Float & err): value(v), l(v-err), h(v+err) {}
    inline EFloat (const Float & v, const Float & l, const Float & h):
        value(v), l(l), h(h) {}
    inline EFloat operator + (const EFloat & b) const {
        return EFloat(value+b.value, prevFloat(l+b.l), nextFloat(h+b.h));
    }
    inline EFloat operator - (const EFloat & b) const {
        return EFloat(value-b.value, prevFloat(l-b.h), nextFloat(h-b.l));
    }
    inline EFloat operator - () const {return EFloat(-value, -h, -l);}
    inline explicit operator float () const {return value;}
    inline Float lower () const {return l;}
    inline Float upper () const {return h;}
    inline EFloat operator * (const EFloat & v) const {
        Float a = l*v.l, b = l*v.h, c = h*v.l, d = h*v.h;
        return EFloat (
            value*v.value,
            prevFloat(min(min(a, b), min(c, d))),
            nextFloat(max(max(a, b), max(c, d)))
        );
    }
    inline EFloat operator * (const Float & v) const {
        return (*this)*(EFloat)v;
    }
    inline EFloat operator / (const EFloat & v) const {
        if(v.l < 0 && v.h > 0) return EFloat(value/v.value, -Infinity, Infinity);
        Float a = l/v.l, b = l/v.h, c = h/v.l, d = h/v.h;
        EFloat ret = EFloat (
            value/v.value,
            prevFloat(min(min(a, b), min(c, d))),
            nextFloat(max(max(a, b), max(c, d)))
        );
        assert(!isNan(*this));
        return ret;
    }
    inline EFloat operator / (const Float & v) const {
        assert(v != 0);
        return (*this)/(EFloat)v;
    }
    inline EFloat & operator += (const EFloat & b) { value += b.value; l = prevFloat(l+b.l), h = nextFloat(h+b.h); return *this; }
    inline EFloat & operator -= (const EFloat & b) { value -= b.value; l = prevFloat(l-b.h), h = nextFloat(h-b.l); return *this; }
    inline EFloat & operator *= (const EFloat & v) {
        Float a = l*v.l, b = l*v.h, c = h*v.l, d = h*v.h;
        value *= v.value;
        l = prevFloat(min(min(a, b), min(c, d))),
        h = nextFloat(max(max(a, b), max(c, d)));
        return *this;
    }
    inline EFloat & operator /= (const EFloat & v) {
        if(v.l < 0 && v.h > 0) {
            value /= v.value;
            l = -Infinity, h = Infinity;
            return *this;
        }
        Float a = l/v.l, b = l/v.h, c = h/v.l, d = h/v.h;
        value /= v.value;
        l = prevFloat(min(min(a, b), min(c, d))),
        h = nextFloat(max(max(a, b), max(c, d)));
        assert(!isNan(*this));
        return *this;
    }

    inline bool operator < (const EFloat & b) const 
    {
        return value < b.value;
    }
    inline bool operator < (const Float & b) const 
    {
        return value < b;
    }
};

bool isNan (EFloat a) {
    return isNan(a.l) || isNan(a.h) || isNan(a.value) || a.l > a.value || a.value > a.h;
}

inline EFloat sqrt(EFloat a) {
    Float v = std::sqrt(a.value);
    Float l = a.l != 0 ? prevFloat(std::sqrt(a.l)) : 0;
    Float h = nextFloat(std::sqrt(a.h));
    assert(!isNan(v) && !isNan(l) && !isNan(h));
    return EFloat(v, l, h);
}

inline EFloat abs (EFloat a) {
    if(a.l >= 0) return a;
    if(a.h <= 0) return -a;
    return EFloat(abs(a.value), 0, a.h);
}

inline bool quadratic (EFloat a, EFloat b, EFloat c, EFloat * x1, EFloat * x2) {
    double de = (double)b.value*b.value-(double)a.value*c.value*4;
    if(de < 0) return false;
    de = std::sqrt(de);
    EFloat s(de, de*float_epsilon);
    EFloat q;
    if ((Float)b < 0)
        q = (b-s)*-0.5;
    else q = (b+s)*-0.5;
    *x1 = q/a;
    *x2 = c/q;
    if(x1->value > x2->value) std::swap(*x1, *x2);
    return true;
}

}

#endif